function vega = Vega(sigma, r, K, S0, dx, dt, T, scheme, type, XminFactor, ...
	XmaxFactor, dsigma)

[plow, vtl, Phi, Delta, Gamma, S] = PDEFD(sigma - dsigma, r, K, S0, dx, dt, T, ...
   scheme, type, XminFactor, XmaxFactor);

[phigh, vtu, Phi, Delta, Gamma, S] = PDEFD(sigma + dsigma, r, K, S0, dx, dt, T, ...
   scheme, type, XminFactor, XmaxFactor);

vega = zeros(length(vtu(:, 1)), 1);

vega(:) = (vtu(:, 1) - vtl(:, 1)) / (2 * dsigma);
 
